Compute Jacobian matrix of Deep learning model?
In the paper The Limitations of Deep Learning in Adversarial Settings, Papernot et. al., 2016, IEEE the chain rule is used, "To express $\nabla F(X∗)$ in terms of $X$ and constant values only."

Earlier is stated, "Our understanding of how changes made to inputs affect a DNN’s output stems from the evaluation of the forward derivative: a matrix we introduce and define as the Jacobian of the function learned by the DNN. The forward derivative is used to construct adversarial saliency maps indicating input features to include in perturbation $\partial X$ in order to produce adversarial samples inducing a certain behavior from the DNN;"

And later, "We define the forward derivative as the Jacobian matrix of the function F learned by the neural network during training. For this example, the output of F is one dimensional, the matrix is therefore reduced to a vector (below). Both components of this vector are computable using the adversary’s knowledge, and later we show how to compute this term efficiently."

$\nabla F(X) = \Big[ \frac {\partial F(X)} {\partial X_1}, \frac {\partial F(X)} {\partial X_2} \Big] \quad\quad (2)$

Between equation (2) and equation (6) your questeion is answered, resulting in an equation where the chain rule has already been applied.

$\frac {\partial F_j (X)} {\partial x_i} = \Big(W_{n+1,j} \, . \, \frac {\partial H_n} {\partial x_i} \Big) \times \frac {\partial f_{n+1,j}} {\partial x_i} \;(W_{n+1,j} \, . \, H_n + b_{n+1,j}) \quad\quad(6)$

It is (6) that you must implement in code, but that is done after the DNN is converged as the paper states above.

Applying the Chain Rule

If $z = f(y)$ and $y = g(x)$, then $\frac {dz} {dx} = f'(g(x)) g'(x)$

Applied twice, if $w = e(z)$, $z = f(y)$, and $y = g(x)$, then $\frac {dw} {dx} = e'(f(g(x)) f'(g(x)) g'(x)$

Applied indefinitely we multiply all the first order derivatives of intermediate values, which has the following first order predicate logic expression.

$\forall \; 0 \le n \lt N \land v_{n + 1} = f_{n + 1}(v_n) \land F_{n+1}(x) = f_{n+1}(F_n(x)) \land g(a) = \dfrac {d} {d a} \, f(a)$
$\exists \; g_n(v_0) = \prod_{\,n = 1}^{\,n &lt; N} {g_n} (F_n(v_0))$

Derivatives of convolution kernel values are the values themselves because they represent attenuation of the form $y = ax + b$.  This is true of first degree linear activation functions too.  The derivative normally used for ReLU (although the x = 0 case is actually undefined) is as follows:

$\dfrac {d} {dx} R(x) = {\begin{cases}x &lt; 0 &amp; 0\\x = 0 &amp; 1\\x &gt; 0 &amp; 1 \end{cases}}$

The derivative of a max function used in max pool (although the u = v case is actually undefined) is as follows:

$\dfrac {d} {dt} M(u, v) = {\begin{cases}u &lt; v &amp; \frac {dv} {dt}\\u = v &amp; \frac {dv} {dt}\\u &gt; v &amp; \frac {du} {dt}\end{cases}}$