text: 4.5. Weight Decay¶
  Now that we have characterized the problem of overfitting and motivated the need for capacity control, we can begin discussing some of the popular techniques used to these ends in practice. Recall that we can always mitigate overfitting by going out and collecting more training data, that can be costly and time consuming, typically making it impossible in the short run. For now, let’s assume that we have already obtained as much high-quality data as our resources permit and focus on techniques aimed at limiting the capacity of the function classes under consideration.
  In our toy example, we saw that we could control the complexity of a polynomial by adjusting its degree. However, most of machine learning does not consist of polynomial curve fitting. And moreover, even when we focus on polynomial regression, when we deal with high-dimensional data, manipulating model capacity by tweaking the degree \(d\) is problematic. To see why, note that for multivariate data we must generalize the concept of polynomials to include monomials, which are simply products of powers of variables. For example, \(x_1^2 x_2\), and \(x_3 x_5^2\) are both monomials of degree \(3\). The number of such terms with a given degree \(d\) blows up as a function of the degree \(d\).
  Concretely, for vectors of dimensionality \(D\), the number of monomials of a given degree \(d\) is \({D -1 + d} \choose {D-1}\). Hence, a small change in degree, even from say \(1\) to \(2\) or \(2\) to \(3\) would entail a massive blowup in the complexity of our model. Thus, tweaking the degree is too blunt a hammer. Instead, we need a more fine-grained tool for adjusting function complexity.
  4.5.1. Squared Norm Regularization¶
  Weight decay (commonly called L2 regularization), might be the most widely-used technique for regularizing parametric machine learning models. The basic intuition behind weight decay is the notion that among all functions \(f\), the function \(f = 0\) is the simplest. Intuitively, we can then measure functions by their proximity to zero. But how precisely should we measure the distance between a function and zero? There is no single right answer. In fact, entire branches of mathematics, e.g., in functional analysis and the theory of Banach spaces are devoted to answering this issue.
  For our present purposes, a very simple interpretation will suffice: We will consider a linear function \(f(\mathbf{x}) = \mathbf{w}^\top \mathbf{x}\) to be simple if its weight vector is small. We can measure this via \(||mathbf{w}||^2\). One way of keeping the weight vector small is to add its norm as a penalty term to the problem of minimizing the loss. Thus we replace our original objective, minimize the prediction error on the training labels, with new objective, minimize the sum of the prediction error and the penalty term. Now, if the weight vector becomes too large, our learning algorithm will find more profit in minimizing the norm \(|| \mathbf{w} ||^2\) versus minimizing the training error. That is exactly what we want. To illustrate things in code, let’s revive our previous example from Section 3.1 for linear regression. There, our loss was given by
  Recall that \(\mathbf{x}^{(i)}\) are the observations, \(y^{(i)}\) are labels, and \((\mathbf{w}, b)\) are the weight and bias parameters respectively. To arrive at a new loss function that penalizes the size of the weight vector, we need to add \(||mathbf{w}||^2\), but how much should we add? To address this, we need to add a new hyperparameter, that we will call the regularization constant and denote by \(\lambda\):
  This non-negative parameter \(\lambda \geq 0\) governs the amount of regularization. For \(\lambda = 0\), we recover our original loss function, whereas for \(\lambda > 0\) we ensure that \(\mathbf{w}\) cannot grow too large. The astute reader might wonder why we are squaring the norm of the weight vector. We do this for two reasons. First, we do it for computational convenience. By squaring the L2 norm, we remove the square root, leaving the sum of squares of each component of the weight vector. This is convenient because it is easy to compute derivatives of a sum of terms (the sum of derivatives equals the derivative of the sum).
  Moreover, you might ask, why the L2 norm in the first place and not the L1 norm, or some other distance function. In fact, several other choices are valid and are popular throughout statistics. While L2-regularized linear models constitute the classic ridge regression algorithm L1-regularized linear regression is a similarly fundamental model in statistics popularly known as lasso regression.
  One mathematical reason for working with the L2 norm and not some other norm, is that it penalizes large components of the weight vector much more than it penalizes small ones. This encourages our learning algorithm to discover models which distribute their weight across a larger number of features, which might make them more robust in practice since they do not depend precariously on a single feature. The stochastic gradient descent updates for L2-regularized regression are as follows:
  As before, we update \(\mathbf{w}\) based on the amount by which our estimate differs from the observation. However, we also shrink the size of \(\mathbf{w}\) towards \(0\). That is why the method is sometimes called “weight decay”: because the penalty term literally causes our optimization algorithm to decay the magnitude of the weight at each step of training. This is more convenient than having to pick the number of parameters as we did for polynomials. In particular, we now have a continuous mechanism for adjusting the complexity of \(f\). Small values of \(\lambda\) correspond to unconstrained \(\mathbf{w}\), whereas large values of \(\lambda\) constrain \(\mathbf{w}\) considerably. Since we do not want to have large bias terms either, we often add \(b^2\) as a penalty, too.
  4.5.2. High-Dimensional Linear Regression¶
  For high-dimensional regression it is difficult to pick the ‘right’ dimensions to omit. Weight-decay regularization is a much more convenient alternative. We will illustrate this below. First, we will generate some synthetic data as before
  representing our label as a linear function of our inputs, corrupted by Gaussian noise with zero mean and variance 0.01. To observe the effects of overfitting more easily, we can make our problem high-dimensional, setting the data dimension to \(d = 200\) and working with a relatively small number of training examples—here we will set the sample size to 20:
  %matplotlib inline import d2l from mxnet import autograd, gluon, init, np, npx from mxnet.gluon import nn npx.set_np() n_train, n_test, num_inputs, batch_size = 20, 100, 200, 1 true_w, true_b = np.ones((num_inputs, 1)) * 0.01, 0.05 train_data = d2l.synthetic_data(true_w, true_b, n_train) train_iter = d2l.load_array(train_data, batch_size) test_data = d2l.synthetic_data(true_w, true_b, n_test) test_iter = d2l.load_array(test_data, batch_size, is_train=False)
  4.5.3. Implementation from Scratch¶
  Next, we will show how to implement weight decay from scratch. All we have to do here is to add the squared \(\ell_2\) penalty as an additional loss term added to the original target function. The squared norm penalty derives its name from the fact that we are adding the second power \(\sum_i w_i^2\). The \(\ell_2\) is just one among an infinite class of norms call p-norms, many of which you might encounter in the future. In general, for some number \(p\), the \(\ell_p\) norm is defined as
  4.5.3.1. Initializing Model Parameters¶
  First, we will define a function to randomly initialize our model
  parameters and run
  attach_grad on each to allocate memory for the
  gradients we will calculate.
  def init_params(): w = np.random.normal(scale=1, size=(num_inputs, 1)) b = np.zeros(1) w.attach_grad() b.attach_grad() return [w, b]
  4.5.3.2. Defining \(\ell_2\) Norm Penalty¶
  Perhaps the most convenient way to implement this penalty is to square all terms in place and sum them up. We divide by \(2\) by convention (when we take the derivative of a quadratic function, the \(2\) and \(1/2\) cancel out, ensuring that the expression for the update looks nice and simple).
  def l2_penalty(w): return (w**2).sum() / 2
  4.5.3.3. Defining the Train and Test Functions¶
  The following code defines how to train and test the model separately on
  the training dataset and the test dataset. Unlike the previous sections,
  here, the \(\ell_2\) norm penalty term is added when calculating the
  final loss function. The linear network and the squared loss have not
  changed since the previous chapter, so we will just import them via
  d2l.linreg and
  d2l.squared_loss to reduce clutter.
  def train(lambd): w, b = init_params() net, loss = lambda X: d2l.linreg(X, w, b), d2l.squared_loss num_epochs, lr = 100, 0.003 animator = d2l.Animator(xlabel='epochs', ylabel='loss', yscale='log', xlim=[1, num_epochs], legend=['train', 'test']) for epoch in range(1, num_epochs + 1): for X, y in train_iter: with autograd.record(): # The L2 norm penalty term has been added l = loss(net(X), y) + lambd * l2_penalty(w) l.backward() d2l.sgd([w, b], lr, batch_size) if epoch % 5 == 0: animator.add(epoch+1, (d2l.evaluate_loss(net, train_iter, loss), d2l.evaluate_loss(net, test_iter, loss))) print('l1 norm of w:', np.abs(w).sum())
  4.5.
  train(lambd=0)
  l1 norm of w: 152.896
  4.5.
  train(lambd=3)
  l1 norm of w: 0.30611762
  4.5.4. Concise Implementation¶
  Because weight decay is ubiquitous in neural network optimization, Gluon makes it especially convenient, integrating weight decay into the optimization algorithm itself for easy use in combination with any loss function. Moreover, this integration serves a computational benefit, allowing implementation tricks to add weight decay to the algorithm, without any additional computational overhead. Since the weight decay portion of the update depends only on the current value of each parameter, and the optimizer must to touch each parameter once anyway.
  In the following code, we specify the weight decay hyperparameter
  directly through the
  wd parameter when instantiating our
  Trainer. By default, Gluon decays both weights and biases
  simultaneously. Note that we can have different optimizers for
  different sets of parameters. For instance, we can have one
  Trainer
  with weight decay for the weights \(\mathbf{w}\) and another without
  weight decay to take care of the bias \(b\).
  def train_gluon(wd): net = nn.Sequential() net.add(nn.Dense(1)) net.initialize(init.Normal(sigma=1)) loss = gluon.loss.L2Loss() num_epochs, lr = 100, 0.003 #}) animator = d2l.Animator(xlabel='epochs', ylabel='loss', yscale='log', xlim=[1, num_epochs], legend=['train', 'test']) for epoch in range(1, num_epochs+1): for X, y in train_iter: with autograd.record(): l = loss(net(X), y) l.backward() # Call the step function on each of the two Trainer instances to # update the weight and bias separately trainer_w.step(batch_size) trainer_b.step(batch_size) if epoch % 5 == 0: animator.add(epoch+1, (d2l.evaluate_loss(net, train_iter, loss), d2l.evaluate_loss(net, test_iter, loss))) print('L1 norm of w:', np.abs(net[0].weight.data()).sum())
  The plots look just the same as when we implemented weight decay from scratch but they run a bit faster and are easier to implement, a benefit that will become more pronounced for large problems.
  train_gluon(0)
  L1 norm of w: 163.57814
  train_gluon(3)
  L1 norm of w: 0.3003651
  So far, we only touched upon one notion of what constitutes a simple linear function. For nonlinear functions, what constitutes simplicity can be a far more complex question. For instance, there exist Reproducing Kernel Hilbert Spaces (RKHS) which allow one to use many of the tools introduced for linear functions in a nonlinear context. Unfortunately, RKHS-based algorithms do not always scale well to massive amounts of data. For the purposes of this book, we limit ourselves to simply summing over the weights for different layers, e.g., via \(\sum_l \|\mathbf{w}_l\|^2\), which is equivalent to weight decay applied to all layers.
  4.5.
  4.5.6. Exercises \mid x) \propto P(x \mid w) P(w)\). How can you identify \(P(w)\) with regularization?