<jupyter_start><jupyter_text>DropoutDropout is a fundamentally new kind of regularization, if you can even call it that. While "the exact mechanism [by which it helps] is unclear," dropout seems to cause individual neurons to be more intrinsically interesting, rather than just co-adapting with their neighbors.The algorithm itself is easy to describe. **During training**, for each mini-batch (or sometimes, each individual example), we designate a random subset of the inputs to each level and simply "turn them off" -- that is, set them to zero.  That means each neuron has to develop without knowing which neurons will actually feed into them.**During testing** (or whenever the neural network is used for actual prediction), we do not drop any inputs off.<jupyter_code>#load_ext autoreload
#autoreload 2<jupyter_output><empty_output><jupyter_text>The CodeRecall our forward propagation code looked like this:<jupyter_code>def forward_prop(weights, biases, activations, input_data):
    l = len(weights)
    
    x = [0] * l # input to level i
    z = [0] * l # un-activated output of level i
    y = [0] * l # activated output of level i
    
    x[0] = input_data
    
    for i in range(0, l):
        expanded_bias = np.ones((x[i].shape[0], 1)) * biases[i]
        z[i] = np.dot(x[i], weights[i]) + expanded_bias
        y[i] = activations[i](z[i])
        
        if i < l-1:
            x[i+1] = y[i]
    
    return x, z, y<jupyter_output><empty_output><jupyter_text>Recall that `x[i]` is the input to level `i`, so `x[0]` is the input to the network, `x[1]` is the input to the first hidden layer, and so on.  We can do all our dropout by modifying `x`:<jupyter_code>def forward_prop_dropout(weights, biases, activations, input_data, drop_rates):
    
    l = len(weights)
    
    x = [0] * l # input to level i
    z = [0] * l # un-activated output of level i
    y = [0] * l # activated output of level i
    mask = [0] * l # dropout mask for the input to level i
    
    # Apply a dropout mask and scaling to x[0]
    mask[0] = np.random.random(input_data.shape) > drop_rates[0]
    x[0] = input_data * (mask[0]) / (1-drop_rates[0])
    
    for i in range(0, l):
        expanded_bias = np.ones((x[i].shape[0], 1)) * biases[i]
        z[i] = np.dot(x[i], weights[i]) + expanded_bias
        y[i] = activations[i](z[i])
        
        if i < l-1:
            mask[i+1] = np.random.random(y[i].shape) > drop_rates[i+1]
            x[i+1] = y[i] * (mask[i+1]) / (1-drop_rates[i+1])
    
    return x, z, y, mask<jupyter_output><empty_output><jupyter_text>We only had to change a few lines -- `8` to keep track of the masks, `11` and `12` to mask the inputs to the network, and `20` and `21` to mask the inputs to the hidden layers. Let's explain the new code:1. `np.random.random(mat.shape)` gives a random matrix matching the given shape. The numbers are uniformly generated from zero to one.2. `mat > p` gives a matrix of the same shape as `mat`, which is True or False, depending on whether the specific entry of the matrix is less than `p`. Underneath, True is 1 and False is 0, which lets us use it for arithmetic.  Note that this sets an entry of matrix to zero with probability `p`, and sets it to one with probability `1-p`.3. `input_data * (np.random.random(input_data.shape) > dropout_rates[0])` multiplies `input_data` by a matrix of zeros and ones, meaning it sets the failures to zero and leaves the rest untouched.4. `mat / (1-p)` increases the values of `mat` by dividing them by `1-p`.  This compensates exactly for the average amount of total input we're losing because of the dropout.  For example if `p` were `1/3`, then we're losing `1/3` of the inputs, so keeping only `2/3` of it.  To compensate, we scale the remaining inputs by `3/2`, so that on average, the total input is the same.  Notes about our ImplementationThere are several small notes worth making in the above.First, by using a completely random matrix in step 1, we are dropping out different inputs in each row.  This means that instead of having a fixed dropout choice for the entire mini-batch, we are dropping out different examples independently.  Although the original paper was ambiguous on this point, further research, such as the the introduction on DropConnect, suggests that this gives more regularization and is better.  It's also convenient to code.Second, we compensate for the dropout by scaling as we dropout.  In the original paper, they did not compensate during the training, but did compensate during testing.  Our method is mathematically equivalent, but has the advantage that the network can be tested (or used) without reference to training parameters like the dropout rate.  What About Training?It doesn't seem fair to penalize weights leading to a neuron which was turned off, and in fact we don't.  This is why we needed to keep track of the masks. The original paper declares "any training example which does not use a parameter will not contribute to the gradient of the parameter."  What this means is that if the neuron's output was disabled in a particular row, disable all contribution of that row to the neuron's weights.Since `back_prop` aggregates the total change for each neuron, but we need to decide example-by-example whether to zero out any changes, we'll have to make the changes within `back_prop`.  The original method looked like this:<jupyter_code>def back_prop(weights, biases, acts, cost_function,
              train_X, train_Y,
              x, y, z):
    L = len(weights) # number of layers
    
    cost_diff = cost_function(y[-1], train_Y, diff=True)
    
    # Gradient of cost at each level
    bp_grad = [0] * L

    # The last level is special
    bp_grad[L-1] = cost_diff * acts[L-1](z[L-1], y[L-1], diff=True)
    
    # The rest of the levels are just gotten by propagating backward
    for i in range(L-2, -1, -1):
        scaled_grad = bp_grad[i+1] * acts[i+1](z[i+1], y[i+1], diff=True)
        bp_grad[i] = np.dot(scaled_grad, weights[i+1].T)
    
    # Now adjust for the weights and biases themselves
    bp_grad_w = [0] * L
    bp_grad_b = [0] * L
    
    for i in range(0, L):
        scaled_grad = bp_grad[i] * acts[i](z[i], y[i], diff=True)

        bp_grad_w[i] = np.dot(x[i].T, scaled_grad)
        
        relevant_ones = np.ones((1, x[i].shape[0]))
        bp_grad_b[i] = np.dot(relevant_ones, scaled_grad)
        
    return bp_grad_w, bp_grad_b<jupyter_output><empty_output><jupyter_text>We can update it with one line of code (not counting passing `masks` as a parameter):<jupyter_code>def back_prop_dropout(weights, biases, acts, cost_function,
                           train_X, train_Y,
                           x, y, z, masks):
    L = len(weights) # number of layers
    
    cost_diff = cost_function(y[-1], train_Y, diff=True)
    
    # Gradient of cost at each level
    bp_grad = [0] * L
    
    # The last level is special
    bp_grad[L-1] = cost_diff * acts[L-1](z[L-1], y[L-1], diff=True)
    
    # The rest of the levels are just gotten by propagating backward
    for i in range(L-2, -1, -1):
        scaled_grad = bp_grad[i+1] * acts[i+1](z[i+1], y[i+1], diff=True)
        bp_grad[i] = np.dot(scaled_grad, weights[i+1].T)
    
    # Now adjust for the weights and biases themselves
    bp_grad_w = [0] * L
    bp_grad_b = [0] * L
    
    for i in range(0, L):
        if i < L-1:
            scaled_grad = bp_grad[i] * acts[i](z[i], y[i], diff=True) * masks[i+1]
        else:
            scaled_grad = bp_grad[i] * acts[i](z[i], y[i], diff=True)

        bp_grad_w[i] = np.dot(x[i].T, scaled_grad)
        
        relevant_ones = np.ones((1, x[i].shape[0]))
        bp_grad_b[i] = np.dot(relevant_ones, scaled_grad)
        
    return bp_grad_w, bp_grad_b<jupyter_output><empty_output><jupyter_text>In the original line `24`, we scale the contributions of the gradient to fit the activiation functions.  Now we also zero it out if that neuron's output was missing; note that we don't dropout from the output, so we need a conditional there.  This change gets propagated across the entire neuron as appropriate in the original lines `26` and `29` (now lines `29` and `32`), as usual.  Testing it OutIn the original paper, they suggest a 20% dropout rate for input units and a 50% dropout rate for hidden units.  They also suggest scaling up the number of hidden units proportionally to compensate for the dropout, as well as dramatically increasing the learning rate and momentum.  Let's try it out:<jupyter_code>from mnist_import import get_mnist_nice

train_X, train_Y, test_X, test_Y = get_mnist_nice()

n = train_X.shape[1]
k = train_Y.shape[1]<jupyter_output><empty_output><jupyter_text>Here's our "optimize" code, which just means spawn a neural network and train it up with a mass of hyperparameters. Note that we now use dropout as an additional hyperparameter:<jupyter_code>def optimize(act_fn, cost_fn, init_fn, learning_rate,
             train_X, train_Y,
             neuron_sizes, num_epochs, batch_size,
             l1_cost=0, l2_cost=0, dropout_rates=None,
             momentum=0
            ):
    np.random.seed(313) # for determinism
    
    # Step 2: initialize
    weights, biases = init_fn(n, k, neuron_sizes)
    acts = [act_fn for _ in range(0, len(weights))]
    acts[-1] = act_sigmoid # last one is always sigmoid
    
    if dropout_rates is None: # None means no dropout
        dropout_rates = [0 for _ in range(0, len(weights))]
    
    # Step 3: train
    t1 = time.time()
    
    weight_velocities = [0 for _ in range(0, len(weights))]
    biases_velocities = [0 for _ in range(0, len(biases))]

    for epoch in range(0, num_epochs):
        # we'll keep track of the cost as we go
        total_cost = 0
        num_batches = 0

        for X_mb, Y_mb in get_mini_batches(batch_size, train_X, train_Y):
            x, z, y, masks = forward_prop_dropout(weights, biases, acts, X_mb, dropout_rates)

            bp_grad_w, bp_grad_b = back_prop_dropout(weights, biases, acts, cost_fn, X_mb, Y_mb, x, y, z, masks)
            l1_grad_w = lasso_cost(l1_cost, weights, biases, diff=True)
            l2_grad_w = ridge_cost(l2_cost, weights, biases, diff=True)

            for i in range(0, len(weights)):
                weight_grad = bp_grad_w[i] / len(X_mb)
                weight_grad += l1_grad_w[i]
                weight_grad += l2_grad_w[i]
                
                weight_velocities[i] = weight_velocities[i] * momentum + weight_grad
                weights[i] -= weight_velocities[i] * learning_rate
                
                biases_grad = bp_grad_b[i] / len(X_mb)
                
                biases_velocities[i] = biases_velocities[i] * momentum + biases_grad
                biases[i] -= biases_velocities[i] * learning_rate

            total_cost += cost_fn(y[-1], Y_mb, aggregate=True)
            num_batches += 1

        cost = total_cost / num_batches # average cost
        print("After {0} epochs, have average cost {1:0.7f}. Took {2:0.3f} seconds so far.".format(epoch, cost, time.time()-t1))
        
        av = np.mean([np.mean(sum(w*w)) for w in weights])
        print("Average norm of each weight is {0:0.3f}".format(av))
        
        _, _, y = forward_prop(weights, biases, acts, train_X)
        train_Y_hat = y[-1]

        train_success = classification_success_rate(train_Y_hat, train_Y)
        print("Got {0:0.3f}% success rate on the training data.".format(100 * train_success))
        
        print()
    
    return weights, biases, acts<jupyter_output><empty_output><jupyter_text>And now let's train it with the dropout methods:<jupyter_code>from basic_nn import *
import time
import numpy as np

act_fn = act_tanh
cost_fn = cost_CE
init_fn = initialize_xavier_tanh

learning_rate = 2

neuron_sizes = [1024]
num_epochs = 40
batch_size = 50

l1_cost = 0.5
l2_cost = 0.25

momentum = 0.95

dropout_rates = [0.2, 0.5]

# Fix up the parameters, as this agrees with what happened in the experiment
l1_cost /= len(train_X)
l2_cost /= len(train_X)
learning_rate *= 1 - momentum

# train it up; note this is the whole dataset now
weights, biases, acts = optimize(act_fn, cost_fn, init_fn, learning_rate,
                                 train_X, train_Y,
                                 neuron_sizes, num_epochs, batch_size,
                                 l1_cost=l1_cost, l2_cost=l2_cost, dropout_rates=dropout_rates,
                                 momentum=momentum
                                )

_, _, y = forward_prop(weights, biases, acts, train_X)
train_Y_hat = y[-1]

train_success = classification_success_rate(train_Y_hat, train_Y)
print("Got {0:0.3f}% success rate on the training data.".format(100 * train_success))

_, _, y = forward_prop(weights, biases, acts, test_X)
test_Y_hat = y[-1]

test_success = classification_success_rate(test_Y_hat, test_Y)
print("Got {0:0.3f}% success rate on the test data.".format(100 * test_success))<jupyter_output>After 0 epochs, have average cost 2.0732929. Took 83.812 seconds so far.
Average norm of each weight is 5.609
Got 89.262% success rate on the training data.

After 1 epochs, have average cost 1.1953199. Took 176.944 seconds so far.
Average norm of each weight is 6.301
Got 90.223% success rate on the training data.

After 2 epochs, have average cost 1.1217761. Took 269.073 seconds so far.
Average norm of each weight is 6.719
Got 90.103% success rate on the training data.

After 3 epochs, have average cost 1.0835093. Took 361.983 seconds so far.
Average norm of each weight is 6.970
Got 91.977% success rate on the training data.

After 4 epochs, have average cost 1.0181313. Took 454.245 seconds so far.
Average norm of each weight is 6.830
Got 92.327% success rate on the training data.

After 5 epochs, have average cost 0.9545442. Took 546.972 seconds so far.
Average norm of each weight is 6.661
Got 93.007% success rate on the training data.

After 6 epochs, have average cost 0.9102258. To[...]<jupyter_text>Note that each epoch is taking a painfully long time to run; this because of two factors:1. We are doing a lot of extra calculations, including generating random matrices, then using these masks for multiplication at various points.2. Our neural networks are a lot bigger, to compensate for the dropout rate.Also, the average costs are a lot higher.  This makes sense - dropout is a lot of regularization, and the neurons never know what they're going to have to work with.  It's a much harder problem than mere stochastic gradient descent, with a lot more noise, but it's working.  Slowly.Finally, we notice that overfitting has been completely defeated, but honestly, so has fitting well.  In the original paper, they described using very large learning rates (10-100 times what was originally useful, which for was was `0.5`) and very large momentum (`0.95` to `0.99`), then described another kind of regularization that complemented these things well, called **max-norm** regularization.  Max-Norm RegularizationThe idea is, for each neuron **n**, we insist that the L2-norm is below some constant *c*.  Recall the L2-norm is the square root of the sum of the squares of the elements.  So the regularization here, rather than integrating it with gradient descent, is a simple scaling of the weights.  The net effect is that the learning rate and momentum sling the weights around violently, then the max-norm regularization sticks them to the surface of the appropriate sphere, so that really we're only changing the angle (but again, sometimes violently).Let's try it out.  To force a vector's norm to be exactly `c`, we can multiply by `c / norm`, which will be less than one if and only if `norm > c`. So, to force a vector's norm to be at most `c`, we can multiply it by `min(1, c / norm)`, which does nothing if `norm  c`.<jupyter_code>def max_norm_scale(mat, bound):
    L2 = sum(mat*mat)**0.5 # compute L2 norms of each column; this is now a row vector
    scale = np.minimum(np.ones(L2.shape), bound/L2)
    return mat*scale
def optimize(act_fn, cost_fn, init_fn, learning_rate,
             train_X, train_Y,
             neuron_sizes, num_epochs, batch_size,
             l1_cost=0, l2_cost=0, dropout_rates=None,
             momentum=0, max_norm=float("inf")
            ):
    np.random.seed(313) # for determinism
    
    # Step 2: initialize
    weights, biases = init_fn(n, k, neuron_sizes)
    acts = [act_fn for _ in range(0, len(weights))]
    acts[-1] = act_sigmoid # last one is always sigmoid
    
    if dropout_rates is None: # None means no dropout
        dropout_rates = [0 for _ in range(0, len(weights))]
    
    # Step 3: train
    t1 = time.time()
    
    weight_velocities = [0 for _ in range(0, len(weights))]
    biases_velocities = [0 for _ in range(0, len(biases))]

    for epoch in range(0, num_epochs):
        # we'll keep track of the cost as we go
        total_cost = 0
        num_batches = 0

        for X_mb, Y_mb in get_mini_batches(batch_size, train_X, train_Y):
            x, z, y, masks = forward_prop_dropout(weights, biases, acts, X_mb, dropout_rates)

            bp_grad_w, bp_grad_b = back_prop_dropout(weights, biases, acts, cost_fn, X_mb, Y_mb, x, y, z, masks)
            l1_grad_w = lasso_cost(l1_cost, weights, biases, diff=True)
            l2_grad_w = ridge_cost(l2_cost, weights, biases, diff=True)

            for i in range(0, len(weights)):
                weight_grad = bp_grad_w[i] / len(X_mb)
                weight_grad += l1_grad_w[i]
                weight_grad += l2_grad_w[i]
                
                weight_velocities[i] = weight_velocities[i] * momentum + weight_grad
                weights[i] -= weight_velocities[i] * learning_rate
                weights[i] = max_norm_scale(weights[i], max_norm)
                
                biases_grad = bp_grad_b[i] / len(X_mb)
                
                biases_velocities[i] = biases_velocities[i] * momentum + biases_grad
                biases[i] -= biases_velocities[i] * learning_rate

            total_cost += cost_fn(y[-1], Y_mb, aggregate=True)
            num_batches += 1

        cost = total_cost / num_batches # average cost
        print("After {0} epochs, have average cost {1:0.7f}. Took {2:0.3f} seconds so far.".format(epoch+1, cost, time.time()-t1))
        
        av = np.mean([np.mean(sum(w*w)**0.5) for w in weights])
        print("Average norm of each weight is {0:0.3f}".format(av))
        
        _, _, y = forward_prop(weights, biases, acts, train_X)
        train_Y_hat = y[-1]

        train_success = classification_success_rate(train_Y_hat, train_Y)
        print("Got {0:0.3f}% success rate on the training data.".format(100 * train_success))
        
        print()
    
    return weights, biases, acts<jupyter_output><empty_output><jupyter_text>The only change required is to insert the max_norm scaling into line `42` of the above. Note that we do not scale the bias (there is no need).  Let's see how this works out.  In the below we use much larger learning rate and momentum, and reduce the L2 penalty to zero, since max-norm regularization does essentially the same job.<jupyter_code>from basic_nn import *
import time
import numpy as np

act_fn = act_tanh
cost_fn = cost_CE
init_fn = initialize_xavier_tanh

learning_rate = 2

neuron_sizes = [1024]
num_epochs = 40
batch_size = 50

l1_cost = 0.5
l2_cost = 0.0

momentum = 0.95

dropout_rates = [0.2, 0.5]
max_norm = 5

# Fix up the parameters, as this agrees with what happened in the experiment
l1_cost /= len(train_X)
l2_cost /= len(train_X)
learning_rate *= 1 - momentum

# train it up; note this is the whole dataset now
weights, biases, acts = optimize(act_fn, cost_fn, init_fn, learning_rate,
                                 train_X, train_Y,
                                 neuron_sizes, num_epochs, batch_size,
                                 l1_cost=l1_cost, l2_cost=l2_cost, dropout_rates=dropout_rates,
                                 momentum=momentum, max_norm=max_norm
                                )

_, _, y = forward_prop(weights, biases, acts, train_X)
train_Y_hat = y[-1]

train_success = classification_success_rate(train_Y_hat, train_Y)
print("Got {0:0.3f}% success rate on the training data.".format(100 * train_success))

_, _, y = forward_prop(weights, biases, acts, test_X)
test_Y_hat = y[-1]

test_success = classification_success_rate(test_Y_hat, test_Y)
print("Got {0:0.3f}% success rate on the test data.".format(100 * test_success))<jupyter_output>After 1 epochs, have average cost 2.8655089. Took 108.901 seconds so far.
Average norm of each weight is 1.868
Got 81.245% success rate on the training data.

After 2 epochs, have average cost 1.6802186. Took 231.707 seconds so far.
Average norm of each weight is 1.939
Got 89.707% success rate on the training data.

After 3 epochs, have average cost 1.1510421. Took 352.698 seconds so far.
Average norm of each weight is 1.992
Got 90.953% success rate on the training data.

After 4 epochs, have average cost 1.0763883. Took 472.746 seconds so far.
Average norm of each weight is 2.025
Got 91.678% success rate on the training data.

After 5 epochs, have average cost 1.0256319. Took 596.486 seconds so far.
Average norm of each weight is 2.005
Got 92.485% success rate on the training data.

After 6 epochs, have average cost 0.9776702. Took 714.953 seconds so far.
Average norm of each weight is 2.004
Got 92.427% success rate on the training data.

After 7 epochs, have average cost 0.9382232. T[...]
{'fraction_non_alphanumeric': 0.09087530198847797, 'fraction_numerical': 0.029641330607693736, 'mean_word_length': 3.045292238301071, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '411997', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 7346, 'n_tokens_neox': 6764, 'n_words': 2546, 'path': '12 - Dropout.ipynb', 'repo_name': 'rodya-mirov/shallow-minds', 'size': '21524'}