Calculating the Jacobian of a neural network - Cross Validated most recent 30 from stats.stackexchange.com 2019-07-17T10:48:48Z https://stats.stackexchange.com/feeds/question/246088 http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://stats.stackexchange.com/q/246088 2 Calculating the Jacobian of a neural network generic_user https://stats.stackexchange.com/users/17359 2016-11-15T16:26:44Z 2016-11-15T19:52:27Z <p>I'm trying to calculate confidence intervals for a neural network (rather than prediction intervals). I'm following <a href="https://www.esa.espci.fr/sites/www.esa.espci.fr/IMG/pdf/2000conf.pdf" rel="nofollow noreferrer">this</a> paper, which treats them in the same framework as any parametric (parameter-involving?) nonlinear model (which a neural net basically is). Calculating these CI's involves computing the Jacobian -- the matrix of partial derivatives of $\hat y$ with respect to the parameters.</p> <p>In the one-outcome regression case, a neural net with two hidden layers (and no biases, WLOG) is $$y = \sigma(\sigma(\mathbf{X\alpha)\beta)\gamma} + \epsilon$$ where </p> <ul> <li>$y$ is the $N\times 1$ outcome </li> <li>$\mathbf{X}$ is the $N\times P_x$ data</li> <li>$\alpha$ is the $P_x \times P_\beta$ first set of weights </li> <li>$\beta$ is the $P_\beta \times P_\gamma$ second set of weights </li> <li>$\gamma$ is the $P_\gamma \times 1$ final set of weights, linking to $y$ </li> <li>$\epsilon$ is the mean-zero error </li> <li>$\sigma()$ is the activation function</li> </ul> <p>The chain rule gives me the following partial derivatives of the fitted model: $$\begin{array}{rcl} \displaystyle\frac{\partial\hat y}{\partial \hat\gamma} &amp; = &amp; \sigma(\sigma(\mathbf{X\alpha})\beta) \\ \displaystyle\frac{\partial\hat y}{\partial \hat\beta} &amp; = &amp; \sigma'(\sigma(\mathbf{X\alpha})\beta)\gamma\sigma(\mathbf{X}\alpha) \\ \displaystyle\frac{\partial\hat y}{\partial \hat\alpha} &amp; = &amp; \sigma'(\sigma(\mathbf{X\alpha})\beta)\gamma\sigma'(\mathbf{X}\alpha)\beta \mathbf{X} \end{array}$$</p> <p>Now, the matrix dimensions on the above expressions are not conformable, which suggests that vectorizing the Jacobian calculation is not possible. Is this true? </p> <p>Next, if I need to calculate the partial derivatives parameter-by-parameter, what is the appropriate way to deal with the network structure when evaluating each derivative?</p> <p>For example, $\alpha_1$ links the $\mathbf{X}_1$ to $\sigma(\mathbf{X\alpha})_1$. That node in turn links to several different upper nodes. Let's say that there are $g$ hidden units in the topmost layer. Should the Jacobian column for $\alpha_1$ include $\displaystyle\sum_g \sigma(\mathbf{X\alpha})_1 \beta_g$ (i.e.: row-wise sum)? Or would the combination of the $\beta$'s take some other form than a sum? </p> <p>Likewise, $\hat\alpha_1$ influences $\hat y$ through all of the top layer hidden-units. So $\hat\alpha_1$'s Jacobian column should also include $\sigma'(\sigma(\mathbf{X\alpha})\beta)\gamma$, treated as a matrix, which evaluates to a $N \times P_\gamma \times P_\gamma \times 1 = N \times 1$ matrix?</p> <p>$\hat\alpha_1$'s Jacobian column would thus be $$\frac{\partial\hat y}{\partial \hat\alpha_1} = \sigma'(\sigma(\mathbf{X\alpha})\beta)\gamma\circ\displaystyle\sum_g\sigma'(\mathbf{X}_1\alpha_1)\beta_g \circ\mathbf{X}_1$$</p> <p>Is this correct?<br> Is there a more efficient way to do this? Can a general statement be made about how much better the analytical Jacobian would be than one calculated by numerical approximation, for example in the <code>numDeriv</code> package in R?</p>