Before the midterm, we were talking about cost functions that could be used to optimize neural
nets.
We introduced the notion of a maximal likelihood estimate for a parameterized family of probability distributions:
\begin{eqnarray*}
\vec{\theta}_{ML} &=& argmax_{\vec{\theta}} p_{model}(\mathbb{X}; \vec{\theta})\\
&=& argmax_{\vec{\theta}}\prod_{i=1}^m p_{model}(\vec{x}^{(i)}; \vec{\theta})
\end{eqnarray*}
As products of small numbers are subject to underflow round off errors, we wanted to develop techniques for computing `\vec{\theta}_{ML}` avoid such errors, so we instead tried to maximize the log of the above.
In an examining this, we showed that such a maximization was equivalent to minimizing the Kullback-Leibler divergence of `hat{p}_{data}` and `p_{mo\de\l}`.
This in turn we argued was equivalent to minimizing the cross entropy:
`-E_{vec{x}~hat{p}_{data}}[log p_{mo\de\l}(\vec{x}; vec{theta})].`
or for labeled data:
`-E_{vec{y}~hat{p}_{data}}[log p_{mo\de\l}(y|\vec{x};vec{theta})]`
We showed in the particular case where the model was a linear Gaussian the above quantity is the same as the mean square error:
`1/m sum_i||y^{(i)} - hat{y}^{(i)}||^2`
We begin today by looking the case where the model is a logistic function.
Cross-Entropy Cost Function - Logistic Function Case
Suppose the values of `y` are in `{0,1}` as we have a binary classifier.
Let `vec{theta}` represent the weights coming into a perceptron.
By using the trick of having a dummy input `x_0` which is always 1, we can imagine the perceptron computes some function `g(vec{theta}^Tvec{x})`.
As before, the cross entropy can be calculated as:
`-E_{\vec{y}~\hat{p}_{data}}[\log p_{mo\d\el}(y|\vec{x};\vec{\theta})] = - 1/m \sum_{i=1}^m \log p_{mo\d\el}(y^{(i)}|\vec{x}^{(i)};\vec{\theta}).`
Assume `g` is the logistic function: `g(z) = 1/(1-e^{-z}) = 1/2 + 1/2tanh(z/2)`.
`g(vec{theta}^Tvec{x})` represents the probability of finding `y=1` and `1-g(vec{theta}^Tvec{x})` represents the probability of finding `y=0`.
For the actual data `y` represents probability of the label being 1, and `1-y` represents the probability of the label being 0.
Since we are taking an average anyway, rather than try to compute `\log p_{mo\d\el}(y^{(i)}|\vec{x}^{(i)};\vec{\theta})` directly we can use in the above the easier to calculate:
`E(\log p_{mo\d\el}(y^{(i)}|\vec{x}^{(i)};\vec{\theta})) = y^{(i)} \log g(vec{theta}^Tvec{x}^{(i)}) + (1 - y^{(i)}) \log (1 - g(vec{theta}^Tvec{x}^{(i)}))`.
This gives the cross entropy loss function as:
`- 1/m \sum_{i=1}^m[y^{(i)} \log g(vec{theta}^Tvec{x}^{(i)}) + (1 - y^{(i)}) \log (1 - g(vec{theta}^Tvec{x}^{(i)}))].`
This approximation is also used for other activation functions as well.
Output Units
The choice of cost function is often related to the choice of the units in our output layer.
For example, suppose we were trying to learn a vector of average statistics `vec{s}` a baseball team might get based on data about features `vec{h}` in a hidden layer computed from players that make up the team.
We might choose the output unit to be a linear unit layer `hat{y} = W^{T}vec{h} + vec{b}`.
Suppose, however, we were trying to predict the value of a binary variable `y`, a binary classification problem.
Using MLE, we define a Bernoulli distribution over `y` conditioned on `vec{x}`.
The net needs to predict `P(y=1|vec{x})`. For this number to be a valid probability it must lie in the interval `[0,1]`.
One way to force this to be true if we used a linear unit might be to define:
`P(y=1|vec{x}) = max(0, min(1,vec{w}^T vec{h} + b))`
that is we clip the value of the linear unit to the desired interval.
This has a drawback when using gradient approaches in that the derivative is 0 outside the unit interval.
Instead, a single perceptron unit as the output could be used with a sigmoid activation function which has a well defined derivative.
Motivating the Sigmoid Function, Softmax
Split the top perceptron of our network into two parts `z = vec{w}^T vec{h} + b` and `y=g(z)` for our activation function. The latter step converts `z` to a probability.
Suppose we want to define a probability distribution over `y` using the value for `z`
We can imagine we initially have an unnormalized probability distribution `tilde{P}(y)` that we divide by an appropriate constant to make a probability distribution.
If we assume that the unnormalized log probabilities are linear in `y` and `z`, we have:
\begin{eqnarray*}
\log \tilde{P}(y) &=& yz\\
\tilde{P}(y) &=& exp(yz),\\
P(y) &=& exp(yz)/ \sum_{y'=0}^{1}exp(y'z)\\
P(y) &=& \sigma((2y-1)z)
\end{eqnarray*}
where `sigma` is the sigmoid function.
This motivates using the sigmoid function in the last layer if one has two possibilities.
If we have a `k`-label classifier, we can generalize this argument to say we should a softmax layer:
`soft\m\ax(vec{z})_i = exp(z_i)/(sum_j(exp(z_j)))`
Training a Softmax Layer
Suppose we want to train a softmax layer with `k` labels.
If training input `vec{x^{(j)}` has label `i`, then the vector `vec{y}^{(j)}` we train with should
be a length `k` vector with all components `0` except the `i`th component which should have value 1.
This kind of training vector is called a 1-hot vector.
Minimizing Functions
We want to update our neural nets weights so as to minimize our cost function.
From calculus, we know that we can find local extrema (local minima or local maxima) of functions by computing their derivatives, setting them to 0, and solving for the minimized variables.
Recall a point `vec{x}` is a local minima for `f` if for some `epsilon >0` and all points `vec{z}` within distance `epsilon` of `vec{x}, f(vec{x})\leq f(vec{z})`.
In learning we often settle for local minima.
It is the second derivative determines whether an extrema is a minima or maxima.
In the case, of a single variable twice differentiable function, `f'(x)=0` and `f''(x) geq 0` implies the point `x` is a local minima.
For multivariable function case, if the second derivative, the Hessian, is positive semi-definite (`vec{x}^TM\vec{x} \geq 0` for all `vec{x} in RR^n` ), then we will be at a local minima.
We will call a function which satisfies this test convex.
If it is not practical to compute the exact form of the derivative then to find extrema, we need to search for them.
Newton Methods For Finding Minima
Recall one way, Newton-Raphson aka Newton's Method (1669), to find the zeros of a function `f(x) = y` is to start with an initial guess of `x_0`, compute the tangent line to `(x_0, f(x_0))`,
see the value `x_1` where that intersects the `x`-axis, and use that as a guess for the zero, and repeat if it is not yet close enough.
This gives the update rule: `x_{n+1} = x_n - f(x_n)/(f'(x_n))`.
Under reasonable conditions this will eventually converge on a zero.
If our cost function `C` were 1 dimensional, one could imagine to minimize `C` we find the zeros of the derivative, `C'` using this method. This would involve `C''` and would assume the function was twice differentiable.
Quasi-Newton methods serve as the basis of one class of minimization techniques in the multi-variable setting.
In the multi-variable setting the above observations mean we need to compute Hessians which often need to be tweaked to be made positive-definite.
One algorithm, Levenberg-Marquardt (1944-1963), often used for surface fitting, requires adding an extra `alpha I` term to the Hessian so that the inverse we need for the Newton update is more likely to behave.
Another class of techniques known as conjugate gradient descent algorithms (Fletcher Reeves 1964, Polak Ribiere 1969) try to avoid computing the inverse Hessian by using the Hessian to compute direction lines to multiply against the gradient.
Nowadays, a simpler, closer to Newton approach known as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (four papers, one for each author, all from 1970) has become more popular in the Neural Net setting.
The BFGS Algorithm
The BFGS Algorithm maintains three quantities: `vec{x}_i`, our current best estimate of a solution; `vec{p}_i`, a search direction; and `B_i`, an approximation to the Hessian. We start from an initial guess `vec{x}_0` and an initial approximation of the Hessian `B_0`. To compute the `i+1` values from the `i` values we then do:
Solve `B_i vec{p} = - grad f(vec{x}_i)` for `vec{p}_i` using an algorithm to solve linear equations (in class we did Gaussian elimination earlier in the semester).
Compute a learning rate (aka step size) for this iteration by searching for `alpha_i = argmin_{alpha} f(vec{x}_i + alpha vec{p}_i)`. Note minimizing `f(vec{x}_i + alpha vec{p}_i)` for `alpha > 0` where everything else is fixed reduces to minimizing a one dimensional function, so could use classic 1-D Newton Raphson, in practice, inexact line search is used with Wolfe conditions used to determine if a sufficient good `alpha` has been found.
Set `vec{x}_{i+1} = vec{x}_{i} + vec{s}_i` where `vec{s}_i = alpha_i vec{p}_i`.
Compute `B_{i+1}` by first defining `vec{y}_i = grad f(vec{x}_{i+1}) - grad f(vec{x}_{i})`, and then setting
`B_{i+1} = B_{i} + (vec{y}_i vec{y}_i^T)/(vec{y}_i^T vec{s}_i) + (B_i vec{s}_i vec{s}_i^T B_i^T)/(vec{s}_i^T B_i vec{s}_i)`
You can think of (1) above as computing an analog of `f(x_n)/(f'(x_n))`.
The linear solving portion of the algorithm above makes it at least quadratic and probably cubic in practice to implement step 1.
There is a variant of this algorithm, limited-BFGS (Nocedal 1980), which is used in practice and avoids the computation of this full `n times n` matrix. For limited-BFGS, the update only takes `O(n)`.
Quiz
Which of the following is true?
All mathematical estimators have some bias.
To obtain an MLE for labeled data, we used conditional probabilities.
The cross entropy function developed for the linear Gaussian case is strictly larger than mean square error.
Gradient Descent
Instead of searching for 0's of the derivative, if we just want to minimize `C(x)`, we can do the simpler procedure of starting with some `x_0`, compute `C'(x_i)`, and use an update `x_{i+1} = x_i - epsilon C'(x_i)`. I.e., we move in the direction that makes the function less by a rate proportional to the current derivative.
The generalization of this method to many dimensions is called the method of gradient descent, or in the case of minimization, the method of steepest descent (Debye 1909). When one has a function `f: RR^m -> RR`, this method's update rule looks like:
`vec{x}_{n+1} = vec{x}_n - epsilon grad f(vec{x}_n)`.
Here the gradient, `grad f(vec{x}_n)`, is a vector of partial derivatives with respect to each basis direction and `-grad f(vec{x}_n)` points in the direction of greatest decrease.
`epsilon >0` here is called the learning rate.
The line search method of choosing the learning rate is to evaluate `f(vec{x} - epsilon grad f(vec{x}))` for several values of `epsilon` and choose the one that results in the smallest objective function value.
Hidden Units
Hidden layers of a neural network are layers not directly connected to the output.
For now, we are assuming the units in our layers are perceptrons and are concerned mainly with what activation function to choose.
Since we are going to train our models by using gradient descent to find a local mimima, it often makes sense to use this as one of the selecting criteria.
A perceptron computes `g(vec{w}^T\vec{x} + b)`. Let `f(vec{x}; vec{w}, b) = vec{w}^T\vec{x} + b` so its partial derivative with respect to a `w_i` or `b` is given by the chain rule:
`(del g)/ (del f) \cdot (del f)/(del w_i)` and `(del g)/ (del f) \cdot (del f)/(del b)` .
Types of Activation
identity function g(z) = z - the resulting perceptron is called a linear unit. `(del g)/ (del f)` is 1 and so the derivative of the `(del g)/(del w_i)` and `(del g)/(del b)` will be consistent with `f` and convergence will tend to be fast.
`g(z)= max{0,z}` - the resulting perceptron is called a rectified linear unit (ReLU).
To compute the derivative software packages usual compute either a left-hand or right hand limit. So even though the derivative is not defined at 0, the handed derivatives are and this usually suffice in practice.
Notice `(del g)/ (del f)` is 1 whenever the unit is active, and so the overall derivatives `(del g)/(del w_i)` and `(del g)/(del b)` will tend to remain large whenever the unit is active resulting in fast convergence.
Also, the second derivative is 0 almost everywhere, so the function will tend to behave in a way that will achieve a local minimum, more so than a complex activation function.
`g(z; alpha) = max(0, z) + alpha min(0,z)` - when `alpha = -1`, this is called absolute value rectification. A unit with a small value for `alpha` is called a leaky ReLU. If `alpha` is a learned parameter, then a unit with this kind of activation is called a parametric ReLU.
Units with `g(z; alpha)` have many of the same benefits as ReLU gates as far as their derivatives are concerned but also have a larger domain over which they are non-zero.
Since this domain determines when gradient descent will work they tend to be slightly more robust.
`g(z)=tanh(z)` or the sigmoid activation - these have the virtue of having everywhere defined derivatives, but as there is only a narrow range where these derivatives differ much from 0, units using these activations tend to train more slowly.
Remarks on Initialization
One thing we haven't mentioned yet is what should be the initial choice of weights we should use when we start to do gradient descent.
If we are using ReLUs then it makes sense that we want the biases `vec{b}` not to be `0` and the weights not to be `0` so that the gradient is non zero and we can get going.
For this reason we often initialize the weight to some random float in a range `[-b, b]` where `b` is usually determined based on the number of neurons in the network (we'll talk more about this later).
