Introduction

In this series of posts, I am trying to develop a workflow for using neural networks in data analysis. So far, I have developed R code for fitting a neural network (NN) by gradient descent, used Rcpp to convert that code to C for increased speed, pictured the search paths and used NNs to simulate datasets for comparing algorithms. Now I want to start defining the workflow.

Neural networks are flexible models that are usually over-parameterised and as a consequence the surfaces created by their loss functions have many symmetries and multiple local minima. Searching these complex surfaces by gradient descent is slow and computationally intensive making it important that the algorithm is tuned to give it the best chance of finding a good solution within a reasonable time-frame.

Six aspects of the algorithm that must be specified by the data analyst at the outset are,

• the architecture
  
• the type of activation function
  
• the scaling of the predictors and responses
  
• the starting values for the weights and biases
  
• the step length (learning rate)
• the stopping rule

The first two choices define the model rather than the algorithm and I will leave those choices for another post. The one aspect of the model choice that I will discuss is the centring of the activation function as this affects the efficiency of the algorithm, but leaves the model essentially unchanged.

The other four choices in the list affect the fitting algorithm rather than the model and are inter-related. Rescaling the data affects both the choice of starting values and the choice of the learning rate. When making these choices, we need to balance the desire to find a good model with the computational burden. Running the gradient descent algorithm from 20 different sets of starting values would obviously be better than starting it from 3 different sets, but the gain is unlikely to be worth the effort.

In this post, I will discuss the algorithmic choices in the context of using a NN as a flexible regression model. The post is presented as an investigation; not everything that I tried worked as well as I had hoped, but I think that the failures are themselves quite informative, so I have left them in.

In this post, I cover a lot of ground in a rather superficial way, but at least it is a first step towards developing a general workflow. There is much, much more that can be said about each of these choices, but you have to start somewhere.

Are we talking hyperparameters?

I came to machine learning after a long career in academic statistics, so I bring with me a lot of baggage. In particular, I draw a distinction between a model, such as logistic regression and a fitting algorithm, such as iteratively reweighted least squares. For me, when the model is hierarchical, the higher level model parameters are referred to as hyperparameters. In machine learning, hyperparameters are (usually) tunable features of the algorithm.

There is, in my opinion, a degree of vagueness in machine learning whereby people talk of improvements to the algorithm as if they were improvements to the model. Complaining about this blurring might seem like nit-picking, but it has consequences. For instance, in my next post, I will ask whether cross-validation (CV) measures the performance of the algorithm or the performance of the model? If you use CV, you need to know which it is.

Centring the activation function

The question that I address here is not whether to use one activation function or another. That choice relates to the family of models that is allowed. Rather, I ask how the chosen activation function should be parameterised in order to make the algorithm as efficient as possible.

A NN creates non-linearity and hence its flexibility through its activation function. The popular sigmoid activation function will demonstrate how this works.

x <- seq(-1, 1, 0.01)
bias <- 1
weight <- 2
y <- 1 / (1 + exp(-(bias + weight*x)))
plot(x, y, type="l", lwd=3, ylim=c(0,1), 
     main = "Figure 2: Sigmoid activation when weight=2 and bias=1",
     xlab="x", ylab="Output")

Scaling the training data

Were we to use inappropriate weights and biases, it would be easy to force all inputs to the sigmoid activation function to be outside of the range (-5, 5) and hence to make the outputs almost flat and uninformative. The message is that a neural network that uses the sigmoid activation function needs weights and biases that keep most of the inputs in the range (-5, 5). At the same time, the range of values taken by the input must not become too small, or once again the output will not change.

NNs are easier to fit when the predictors in the input layer are all on similar scales and when that scale matches the outcome of the activation function. That way, weights and biases throughout the network will be applied to quantities of a similar magnitude. What we want to avoid is the situation in a weight in one part of the network is say, 1e-5 and a weight elsewhere is 1e9. This would lead to derivatives of greatly varying size and a numerically unstable algorithm. It would also make it difficult to select a step length/learning rate that is suitable for the whole network.

Since the zero-centred sigmoid has an output range of (-0.5, 0.5), it would make sense to scale the predictors in the same way. The exact range is not critical; anything close to (-0.5, 0.5) would do. A simple transformation that would work for well-behaved data is

\[ \frac{x - \text{mid-point}(x)}{ \text{range}(x)} \]

Unfortunately, this transformation would not be very robust. Remember that we are creating an automatic algorithm, so we do not want to inspect each predictor to ensure that there are no outliers distorting the range. Instead, we might use

\[ \frac{x - (\text{Q}_{95}(x) + \text{Q}_{05}(x)) / 2}{ \text{Q}_{95}(x) - \text{Q}_{05}(x)} \]

where \(Q_{05}\) is the 5% quantile. The range of the transformed predictor will be slightly greater than (-0.5, 0.5), but the vast majority of values will be within the target range.

It is common, at least in regression problems, not to apply an activation function to the final output of the network, that way the output can take whatever value is appropriate to the response of the particular problem. To help make all weights and biases have a similar order of magnitude, the final outputs to the network will need to be scaled so that they too have a value in the range (-5, 5). A possible robust transformation is,

\[ \frac{10 \ (y - (\text{Q}_{95}(x) + \text{Q}_{05}(x)) / 2)}{\text{Q}_{95}(y) - \text{Q}_{05}(y) } \]

Starting Values

The requirement that inputs to a zero-centred sigmoid activation function are in (-5, 5) while the outputs are in (-0.5, 0.5), suggests the use of weights in the range (-10, 10). However, when there are many inputs to a node, the weight given to each one needs to be reduced. Suppose that a particular node has L inputs, we might expect the weights to be in the range (10/L, 10/L). In fact, this is a little conservative because there may be some cancellation with some of the inputs being positive when others are negative. A sensible compromise is, \[ \left( \frac{-10}{\sqrt{L}} , \frac{10}{\sqrt{L}}\right) \]

This type of weight adjustment has been suggested many times based on more mathematical derivations, for example, He et al ( arXiv:1502.01852) suggested that the weights should be reduced by the square root of L for the ReLU activation function and Glorot and Bengio (“Understanding the difficulty of training deep feedforward neural networks.” International conference on artificial intelligence and statistics. 2010) suggested reducing by the root of the number of inputs plus the number of outputs.

The biases are added to the weighted combination of the inputs and the final result needs to be in (-5, 5), so certainly most biases will need to be within this range.

Neural networks are full of symmetries, for instance we could swap any two nodes in a layer provided that we also swap their input and output parameters. By having a good variety of starting weights and biases, you increase the chance of hitting on something close to one of the symmetrical minima. The general rule is that convergence will be quicker when the starting values are not all similar.

Step length (Learning rate)

I am torn between referring to the step length, as one would in statistics, and the learning rate, as one would in machine learning. In the context of gradient descent they mean the same thing, but step length is a little more descriptive, so I will go with that.

Choosing the step length is not trivial, but there are some intuitions that might help guide us

• too small a step length and progress towards the minimum will be slow
  
• too large a step length and the algorithm might jump over the minimum leading to oscillation or even divergence
  
• it makes sense to start with a large step length and to decrease it as we approach the minimum
  
• a good step length for one parameter might not be ideal for other parameters, although scaling the data will greatly reduce this problem
  
• the more parameters in the model the smaller the step length, because each iteration changes all of the parameters simultaneously
  
• a small step length will enable the exploration of small features in the surface of the training loss. However, such narrow features are likely to be specific to the particular training dataset and so might not generalise. Broad minima are more important because they reflect trend rather than noise.

Common-sense intuitions are useful, but they are not evidence, and we need to keep in mind that there will be times when even the best intuitions fail us.

Stopping Rules

With neural networks, the gradient descent algorithm tends to progress through cycles with periods of rapid improvement followed by longer periods when progress almost plateaus. You do not want to stop the algorithm during a plateau, if a period of rapid improvement is still to be reached. At the same time, the algorithm cannot be allowed to go on forever in the vain hope of improvement.

The stopping rule could be based on changes to the training loss, or to changes to the parameters, or to the sizes of the derivatives. When automatic step length adjustment is used, it is important that a minimum step length is specified, so that the algorithm does not make countless negligible changes.

So far the discussion has assumed that the goal is to stop when the algorithm reaches a minimum in the training loss. As I explained in my post on search paths, the real aim is to minimise the expected loss and not the training loss. I will consider early stopping based on expected loss in my next post, which is on cross-validation.

Another limit of this investigation is that I assume that the neural network is not so over-parameterised that the loss can be reduced to zero. Such models can exhibit strange behaviour, such as a test loss that increases and then comes down again. Another topic for a future post.

Possible strategies

In my investigation, I planned to use simulated data to compare two strategies

Comparison 1: The Wave

My first comparison is based on a simple regression. Using the method described in my previous post, I used a (1, 4, 8, 4, 1) NN with different randomly generated weights and biases and a standard sigmoid activation function to create a simulated dataset. Figure 3 shows the relationship between the input and the output and a single set of 100 training values. I call it a wave to make it easier to refer to.

Comparison 2

My second comparison is based on a single dataset of n=1000 observations with 6 inputs (predictors, X) and one output (response, Y). The data were generated in such a way that the correct model would have an expected RMSE loss of 1. The training data are shown in figure 12.

These data are analysed with a (6, 10, 1) neural network. Option A and B are run as originally defined except that I ran 50,000 iterations and a stopping rule for Option B based on a 1e-5 percent change in the loss. Figure 13 shows the patterns of loss by iteration. Option B converges sooner, as we might have anticipated from the first comparison.

The test loss of the fitted models is quite a bit lower with Option B suggesting that Option A has converged to a slightly worse solution for all three sets of starting values. This can be seen more clearly from the search paths in figure 14.

Summary

This has been a rather meandering investigation that has raised as many questions as it has answered. Here are the points that I take from it.

• some scaling of the predictors and responses is essential
  
• starting values need to be ball-park sensible so that inputs stay within the active part of the activation function
  
• the exact choice of starting values and step length affects the speed of convergence more than the final solution
  
• models that look very different can actually be equivalent - there are many symmetries in the loss function of a neural network
  
• Option B performed much better than Option A

I am very aware that these conclusions need to be confirmed in a much greater range of scenarios, but I suspect that Option B would prove to be sensible for a wide range of problems. I will use it in future posts and refine it as necessary.

In practice, the time taken by the gradient descent algorithm means that it is only suitable for small sets of training data; stochastic gradient descent (SGD) works better for larger datasets and so it is not worth putting too much effort into refining the gradient descent workflow, when the real need is for a workflow for SGD.

Appendix: Code Changes

The C code used in this post can be found on my GitHub pages as cnnUpdate02.cpp.

I added options to cfit_nn() and cfit_valid_nn() to give the use the possibility to

• use of the zero-centred sigmoid
  
• report the loss to a specified schedule
  
• use scheduled step length reduction
  
• automatically halve the step length if the loss increases
  
• calculate the test loss to a specified schedule
  
• stop early when changes in the training loss are small

The C function Lin() returns the number of inputs to the node served by each weight.