5.1 Supervised Learning

Accordingly, let us assemble all the samples into two matrices:⁴⁴ 4 These are the transposes of the matrices used in most developments of linear regression, but they are consistent with all the other standard conventions used throughout this book.

In place of the preceding models, let us simply look for a solution to the linear equation

The equation has a unique solution if and only if $\mathbf{{Y}}$ is square—that is, there are precisely as many samples (${N}$) as dimensions (${M}$)—and full-rank. When $\mathbf{{Y}}$ is “tall” (more dimensions than samples, ${M}>{N}$), the solution is underdetermined. We have lately addressed this with regularization. When $\mathbf{{Y}}$ is “fat” (more samples than dimensions, ${N}>{M}$), the solution is overdetermined: no matrix $\mathbf{{D}}$ satisfies Eq. 5.9, and a notion of “best fit” must be imposed in order to select just one of the many approximately satisfactory matrices.

Still, let us proceed somewhat naïvely and simply look for an obvious linear-algebraic solution. If $\mathbf{{Y}}$ is full-rank in addition to being fat, then the Gram matrix (in sample space) $\mathbf{{Y}}\mathbf{{Y}}^{\text{T}}$ is square (${M}\times{M}$) and full-rank, and therefore invertible. Hence, starting from Eq. 5.9,

Again we have recovered the normal equations, this time apparently with even fewer assumptions. Lest the normal equations seem inevitable, let us recall that the choice to right multiply by the matrix $\mathbf{{Y}}^{\text{T}}$ was, although perhaps obvious, also arbitrary. Note in particular that, under our assumption that $\mathbf{{Y}}$ is fat, choosing $\mathbf{{D}}$ according to the normal equations does not actually satisfy Eq. 5.9—indeed, no choice of $\mathbf{{D}}$ could, because the problem is overdetermined. What notion of “best fit” have we tacitly imposed in order to select one value for $\mathbf{{D}}$?

The answer can be read off the derivation in Eq. 5.2 above: The normal equations arise from squared-error penalties, and accordingly the procedure is sometimes called the method of least squares^††margin: method of least squares . Evidently the assumption of normally distributed errors and the least-squares penalty are two sides of the same coin. In the case of a linear map, when $\mathbf{{Y}}$ is full rank, the least-squares solution exists. And since the quadratic is convex, the solution is unique; that is, the normal equations yield the global minimizer of the least-squares penalty. As a “sanity check,” we can also reassure ourselve that when Eq. 5.9 really does hold, i.e. when $\mathbf{{Y}}$ is invertible, the normal equations reduce to the unique solution, $\mathbf{{D}}=\mathbf{{X}}\mathbf{{Y}}^{-1}$.

General form:

Cross-entropy minimizing parameters:

The optimal moment parameters occur when the gradient with respect to the natural parameters is zero….

Generalized linear model, assuming canonical response function:

Gradient:

When $\bm{\mu}_{\check{x}|\check{y}}$ is a linear function of the parameters, $\mathbf{{D}}$, as for Gaussian noise, this gradient can simply be set to zero to solve for the optimal parameters.

When it is not, we need to adopt an iterative approach. Here we will use the Newton-Raphson algorithm discussed above. This requires the Hessian; for simplicity, we write it for scalar outputs only:

Then the IRLS update is

This is the solution to (i.e., the normal equations for) a weighted least-squares problem—cf. Eq. 5.10—but with the parenthetical quantity, rather than ${X}$, as the output variable. To generate some intuition about what precisely this quantity is, let us derive IRLS from a different perspective.

We return to the cross-entropy gradient, Eq. 5.13. The problem with setting this gradient to zero and solving for the optimal parameters is, as noted above, that $\mu_{\check{x}|\check{y}}$ is not generally a linear function of those parameters, $\mathbf{{D}}$. However, if we concede that we must take an iterative approach, we can perhaps approximate this function as a linear displacement from the previous value of the parameters. That is, letting $\lambda$ denote the inverse canonical link function, $\mu^{(i)}_{\check{x}|\check{y}}=\lambda({\bm{d}^{(i)}}^{\text{T}}\bm{y})$, we approximate

Substituting this into the gradient, Eq. 5.13, we have

This matches the update generated with Newton-Raphson, Eq. 5.15. Thus, IRLS can be seen as approximating the gradient as locally linear in the parameters; or, perhaps more intuitively, as approximating the loss (the cross entropy) as locally quadratic in these parameters. Accordingly, we seek a least-squares solution (because of the quadratic cost), but apply it iteratively (because the cost is an approximation). The loss differs at each iteration because the “outputs” change, but also because the data are reweighted by $\sigma^{2}_{\check{x}|\check{y}}$. This is a reweighting, rather than just a weighting, because ${\sigma^{2}_{\check{x}|\check{y}}}^{(i)}$ depends on the current parameters, $\bm{d}^{(i)}$.

Having derived IRLS, we get a feel for the algorithm by applying it to a few cases, beginnning with Poisson noise. For simplicity, we consider a single, scalar “output” or dependent variable; and the canonical link function, the natural logarithm. This makes the mean of the distribution an exponentiated (inverse canonical link) linear function of the “input” or independent variables, ${\bm{Y}}{}$, i.e. $\mathbb{E}_{{\check{X}}{}|{\bm{\check{Y}}}}{\mathopen{}\mathclose{{}\left[{\check{X}}\middle|{\bm{Y}}{}}\right]}=\exp\mathopen{}\mathclose{{}\left\{\bm{d}^{\text{T}}{\bm{Y}}}\right\}$. Recall that for Poisson random variables, the variance is equal to the mean. Inserting these values into the first line of Eq. 5.15, yields

We consider the case where the scale parameter is known, but the shape parameter is not. The reverse is much more commonly considered (it is the “gamma” GLiM built into Matlab’s glmfit, for example), probably because the sufficient statistic for the shape parameter is just the sum over independent output samples, $\sum_{n}\check{x}_{n}$. For simplicity, we consider outputs of length one (${K}=1$) only, which makes the parameters a (row) vector, $\bm{\theta}^{\text{T}}$. Now, the log-likelihood of the parameters under the gamma distribution is:

Here we have made explicit that the only natural parameters depending on the inputs—and the parameters—are those associated with the shape parameter: $\eta^{n}_{k}=k_{n}-1$. The scale parameter is assumed constant across all samples. Since $k_{n}$ must be positive, we define the link function as:

The trailing 1 only clutters the derivation, so we work instead with $k_{n}$, starting with its gradient with respect to the parameters:

We shall also need the first and second the derivatives of the log-partition function with respect to $\eta^{n}_{k}$, i.e., the expected value and variance of $\log\check{x}_{n}$:

with $\psi^{(i)}(\cdot)$ the $i^{\text{th}}$-order polygamma function. Putting these pieces together, the gradient (with respect to the parameters) of the entire cost function above is

where $v_{n}=(\log\check{x}_{n}-\mu_{n})k_{n}$.

The Newton-Raphson algorithm also requires the Hessian. Differentiating a second time we find

with $\mathbf{\Lambda}$ a diagonal matrix with entries $d_{nn}=-k_{n}^{2}\sigma^{2}_{n}+k_{n}\log\check{x}_{n}-k_{n}\mu_{n}$. Notice, however, that if we use the expected Hessian, the term $\log\check{x}_{n}-\mu_{n}$ vanishes, leaving $d_{nn}=-k_{n}^{2}\sigma^{2}_{n}$.

It is another question how to set the parameters $\mathbf{{D}}_{l}$ and $\bm{d}_{l}$ at every layer in order to approximate any particular function. Part of the price of arbitrarily complex functions is non-convexity and expensive Hessian computations, and accordingly we set aside Newton-Raphson and focus on first-order methods. Computing the cross-entropy gradient with respect to all the parameters is conceptually straightforward, but to do so efficiently will require some punctilious bookkeeping. This is considerably simplified by differentiating with respect to matrices and vectors rather than scalars; the relevant rules are derived in Section B.1.

As throughout this chapter, we descend the gradient of the conditional cross-entropy. But to focus on the new aspects of this process introduced by the composition of functions, let us ignore the specific distribution family and refer to the surprisal simply as $\mathcal{L}$:

With some foresight, let us begin with the gradients with respect to the biases. Interpreting $\mathcal{L}$ as a function of $\bm{\check{z}}_{{L}}$, and then applying the chain rule⁵⁵ 5 Recall that the chain rule for derivatives with respect to column vectors accumulates terms on the left, rather than the right; see Section B.1 through Eq. 5.17, we find the gradient at the final layer to be

where $\mathbf{J}_{{L}}$ is the Jacobian of the nonlinearity, $\bm{\psi}_{{L}}$, at layer ${L}$. Because $\bm{\psi}_{{L}}$ acts elementwise on its vector argument, $\mathbf{J}_{{L}}=\mathbf{J}_{{L}}^{\text{T}}$ is diagonal. We have assigned a symbol to this gradient in anticipation of reusing it.

Proceeding to the penultimate layer, we find

and at the third-last layer,

This pattern now repeats at subsequent layers. Hence the sequence of bias gradients can be computed with the backward recursion^††margin: implementation note: Because $\mathbf{J}_{l}$ is diagonal, these matrix-vector products can be computed efficiently as elementwise vector products.

Next, we take derivatives with respect to the weight matrices, beginning with the output layer (cf. Eq. B.7):

Proceeding to the penultimate layer, we find:

Proceeding backwards to earlier layers continues the pattern, so the weight gradients can be computed with the same backward recursion, Eq. 5.18, and the formula

Eqs. 5.18 and 5.19 define the backpropagation algorithm^††margin: backprop for a generic neural network, although applying it to any individual case requires specifying $\bm{\psi}_{l}$ and ${\check{p}\mathopen{}\mathclose{{}\left(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\check{x}}{}\middle|\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{y}{}}\right)}$, which we have left generic here. The weights and biases are updated simply by by moving downhill along the gradients of the cross entropy, which are the expected values of the gradients just derived:

Hence,

although more sophisticated forms of gradient descent are typically used in practice.

Suppose we have reason to believe that some of the structure of our network is “conserved” or repeated, as in Fig. LABEL:fig:XXX. For instance, we might imagine segments of the network to correspond to slices of time, in which case it is reasonable to treat every time slice as identical in structure, even though the inputs and outputs change over time. Or we might have reason to believe that certain spatial structures are repeated. In either case, the concrete meaning of the assumption is that the same parameters show up in multiple places in the network. The upshot for backprop is that a little care must be taken in distinguishing partial from total derivatives, after which the algorithm goes through as in the previous section.

To make this concrete, we consider a recurrent neural network (RNN)^††margin: recurrent neural network (RNN) with inputs ${\bm{Y}}_{1},\ldots,{\bm{Y}}_{{T}}\sim p(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{y}_{1},\ldots,\bm{y}_{{T}})$ and “hidden”-unit activities $\bm{\check{z}}_{t}$ recursively computed as

For simplicity, we let the recurrent layer also be the output layer—i.e., the layer at which the loss function is evaluated—in which case the most general loss function can be expressed as $\mathcal{L}(\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\check{z}}_{0},\ldots,\leavevmode\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\bm{\check{z}}_{T},\bm{\theta})$.⁶⁶ 6 We call this the loss function, but the complete loss is the average of this function under the data. Notice that we let it depend on the hidden units at every time step, as well as directly on the parameters. To compute parameter changes, we take the (total) derivative of the loss function with respect to those parameters:

The first equality is just standard calculus; the second follows because $\bm{\theta}$ only occurs in the final terms in the chains of derivatives.⁷⁷ 7 In a common abuse of notation, I have written partial derivatives with respect to variables, rather than functions—otherwise the arguments would have to be included to disambiguate the derivatives of $\bm{\psi}$ at different moments in time. Thus, one can either compute the total effect of $\bm{\theta}$ on each $\bm{\check{z}}_{t}$, and then compute the direct effects of the latter on $\mathcal{L}$; or one can compute the direct effect of $\bm{\theta}$ on each $\bm{\check{z}}_{t}$, and then compute the total effect of the latter on $\mathcal{L}$.^††margin: Exercise LABEL:ex:partialtotal

Because of Eq. 5.21, the total effect on $\mathcal{L}$ of $\bm{\check{z}}_{t}$, the current hidden state, is its direct effect, plus the effects caused by its influence on future hidden states, $\bm{\check{z}}_{t+k},k>0$. This can be expressed recursively as:

The recursion is initialized at $\mathop{}\!\mathrm{d}{\mathcal{L}}/\mathop{}\!\mathrm{d}{\bm{\check{z}}_{T}^{\text{T}}}=\partial\mathcal{L}/\partial\bm{\check{z}}_{T}^{\text{T}}$ and run backwards in time.

The backprop-through time algorithm consists of Eqs. 5.22 and 5.23. Applying it to a particular model requires working out three partial derivatives, $\partial\bm{\check{z}}_{t+1}/\partial\bm{\check{z}}_{t}^{\text{T}}$, $\partial\mathcal{L}/\partial\bm{\check{z}}_{t}^{\text{T}}$, and $\partial\mathcal{L}/\partial\bm{\theta}^{\text{T}}$, for a particular choice of recurrent function and loss function. Using the standard recurrent function for a neural network (cf. Eq. 5.17):

and applying Eqs. 5.22 and 5.23 yields: